Quasi-stationary evolution of systems driven by particle evaporation 
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Abstract 

We study the quasi-stationary evolution of systems where an energetic confinement is unable to completely retain their 
constituents. It is performed an extensive numerical study of a gas whose dynamics is driven by binary encounters and its 
particles are able to escape from the container when their kinetic energies overcome a given cutoff it c . We use a parametric 
family of differential cross sections in order to modify the effectiveness of this equilibration mechanism. It is verified that 
when the binary encounters favor an effective exploration of all accessible velocities, the quasi-stationary evolution is reached 
when the detailed balance is imposed for all those binary collisions which do not provoke particle evaporation. However, the 
weakening of this effectiveness leads to energy distribution functions which could be very well fitted by using a Michie-King-like 
profile. We perform a theoretical analysis, in the context of Hamiltonian systems driven by a strong chaotic dynamics and 
particle evaporation, in order to take into account the effect of the nonhomogeneous character of the confining potential. 

PACS numbers: 02.50.-r; 05.20.-y 



I. INTRODUCTION 

Gravitation is a generic example of a selfgravitating in- 
teraction which is not able to effectively confine particles: 
it is always possible that some of them reach the suffi- 
cient energy for escaping from a putative system, so that 
the astrophysical systems undergo an evaporation pro- 
cess, and therefore, they will never be in a real thermo- 
dynamic equilibrium. This is precisely the origin of the 
long-range singularity of the thermodynamical descrip- 
tion of astrophysical systems PISHjSIEISHHjU E3 • 

On the other hand, the N-body gravitating systems 
satisfy a strong chaos criterion (see Ref . [Til] ) . supporting 
the assumption of an increasingly uniform spreading of 
orbits over the constant energy surface or its intersection 
with any other integral of motion as the angular momen- 
tum. The competition of these dynamical behaviors pro- 
duces the achievement of a quasi- stationary regime in the 
evolution of the astrophysical systems, allowing in this 
way the applicability of some ordinary equilibrium sta- 
tistical mechanics methods [Illllllllllillllll[3[ll| 

The aim of the present paper is to study the quasi- 
stationary evolution of systems driven by a conservative 
microscopic dynamics which is also strongly affected by 
particle evaporation. Our focus will be mainly on the 
analysis of the effect of evaporation under the modifi- 
cation of the efficiency of the equilibration mechanisms 
(sections ITT1 and ITTT1 ) . The interest in such study is justi- 
fied as follows. 
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As discussed elsewhere, the equilibration mechanisms 
are closely connected with the chaotic properties of the 
microscopic dynamics. This last ones explain both, qual- 
itatively and quantitatively, whether or not the dynam- 
ics provides an effective exploration of all those accessi- 
ble microscopic configurations, and consequently, of how 
strong is the system tendency towards its equilibration 
[Hill 13 13 13 El- Hence, it is clearly evident that 
the influence of a dissipative process, as the particle evap- 
oration, is significantly affected by how efficient are the 
system's equilibration mechanisms. 

On the other hand, as a by-product of the above study, 
we will perform in section I1VI a theoretical analysis of 
the effect of evaporation on the distribution function over 
the single particle six-dimensional /i-phase space, / (r, p). 
This analysis is carried out in the context of Hamiltonian 
systems driven by a very strong chaotic dynamics, and 
our interest is the consideration of the nonhomogeneous 
character of the confining interaction. The result is con- 
trasted with the requirements of the Jeans theorem [T?) 
for the collisionless selfgravitating systems, as well as to 
the Michie-King model of globular clusters and elliptical 
galaxies [3 EI 13 El ■ 



II. A SIMPLE MODEL 

It was done in the past an important set of works de- 
voted to the analysis of particle evaporation in the con- 
text of astrophysical systems in the framework of binary 
encounters (see for example in 0, |3 E3 an d ref- 
erences therein). We are interested in performing a nu- 
merical closer look at of these studies by considering the 
binary collisions from a perspective which allows us to 
control the effectiveness of this equilibration mechanism. 
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Considering the general lines of the classical paper of 
Spitzer and Harm we assume that the following 

idealizations procee: The potential energy is constant 
within the volume of the system with spherical symme- 
try but vanishes in the outside region: 



u (r) 



—u c r < R 
r>R 



(1) 



This first idealization disregards the effects of the non- 
homogeneous character of the confining interaction, and 
therefore, this self-gravitating model is equivalent to a 
gas of particles enclosed in a special container which al- 
lows a given particle to escape when its kinetic energy 
overcomes a given cutoff u c . The second idealization is 
to assume that the particle mean-free-path is many times 
the characteristic size (or length scale) of the system, 
which allows us to take into account only binary colli- 
sions. Consequently, a given particle that gains sufficient 
velocity after a collision will evaporate, since the proba- 
bility that this particle loses the excessive energy before 
escaping from container is very small. The third ideal- 
ization is to suppose that the system is embedded in a 
bidimcnsional space. 

The following toy dynamics will be taken into account 
in order to perform a numerical study of the above model. 
Let us consider a system with TV particles and assign to 
each of them a vector velocity v and a binary number 
6, being and bk the corresponding data for the A;-th 
particle. The binary number 6 = indicates that the 
corresponding particle has escaped from the system, or 
to the contrary that it remains trapped when 6=1. The 
dynamical evolution of the model system is carried out 
as follows: 

• For each computational step, a sample of M inter- 
acting pairs are chosen one by one at random. 

• The pair interaction is suppressed when at least one 
particle has already escaped from the system. 

• When the pair interaction is possible, the particles 
change their velocities as follows: 
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being n an unitary vector defined by: 
n = [cos (0 O + 0) , sin (9 Q + 9)] , 



(2) 



(3) 



where the phase 9 takes a random value belong- 
ing to the interval (— tt, tt), being 9q the phase of 
the relative velocity v = (vi — v 2 ) /2. The above 
equation ensures the energy and linear momentum 
conservation of the pair during the collision. If after 
the collision the kinetic energy of a given particle 
of the pair is greater than the escape energy it c , 
its corresponding binary number changes to zero, 
which means that this particle has escaped from 
the system. 
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FIG. 1: Effect of the deformation parameter a in the angle tp 
in which is enclosed the 50 % of deflections during the binary 
encounters, that is, tp = atan [0.5 arctan (~)1 . Dash line rep- 
resents the asymptotic value for the uniform differential cross 
section ©■ 



• The random phase 9 is distributed according to the 
following generic differential cross section: 
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where a is the total cross section, a, a positive 
real parameter which drives the form of the in- 
teraction, and K, a normalization constant, K = 
arctan (it /a). 

The specific form of the differential cross section Q 
allows us to carry out a comparative study of the quasi- 
stationary evolution features when the efficiency of the 
equilibration mechanism is modified: the variation of the 
deformation parameter a affects the way that the binary 
encounters allow the particles to explore the velocities 
which are stable under the evaporation. 

The FIG n shows the effect of the deformation param- 
eter a on the angle tp in which is enclosed the 50 % of the 
binary encounters, \9\ < tp. When a is large enough, the 
expression @ becomes in a uniform distribution: 



da (9) 



d9_ 



(5) 



which will be used as a reference. 

Let Ni be the number of particles remaining in the 
system. We consider the p = N\/N as the relative den- 
sity of the gas. We denoted by e, the kinetic energy per 
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particle of all those particles remaining in the system, 
that is: 



1 N 1 

1 fc=i 



(6) 



The probability that the interaction takes place for a 
given pair is: 




Pcoll 



Ni {Ni - 1) 
N(N-l) 



P 



(7) 



being Ni and N large enough. The expectation value 
of collisions which take place in a computational step, c, 
can be estimated by c ~ Mp co \\ ~ Mp 2 . 

According to the kinetic theory of gases, the rate of 
collision events per the time unit, uj, can be estimated as 
uj ~ ^an 2 vV, being a the total cross section of the par- 
ticles, n, the particles per unit of volume, v, the expecta- 
tion value of the relative velocity among the particles of 
the system, while V is the volume of the gas container. 
Let At be the temporal increment in each computational 
step. Since c ~ toAt, a direct comparison allows us to set 
the equivalences: p ~ n, N ~ V, and therefore, the num- 
ber of samples is each computational step should grow 
proportional to the system size, that is, M ~ pN , where 
p = ^anovAt, being no the initial density of particles 
per unit of volume. If m represents the particles mass, 
the characteristics units for velocity and time are the fol- 
lowing: v c — y // 2u c /m, To — 2/ariQVc, where v c is also the 
cutoff velocity. The temporal increment in each compu- 
tational step At can be estimated by using the average 
relative velocity among the particles of the system as 
At = Xqp/v, being: 



FIG. 2: Geometrical visualization of the inequality l|15jl . 



of n (v, r) in the next computational step r + 5t can be 
expressed as follows: 



— n (v, t)= p [I in (v; t) - I out (v; r)] , (9) 

being 7j„ (v; r) the probability of some particles to reach 
a final velocity v after a collision, while I ou t (v; r) repre- 
sents the probability of some particles with initial veloc- 
ity v to experience a collision in the next computational 
step. These probabilities are given by the formulae: 



Iin ( v ; r ) = J J d 2 v 1 d 2 v 2 K(v\ vi, v 2 ) n (v 1; t) n(v 2 ,r) , 

SxS 

(10) 



lout (v; r) = n (v, t) / d 2 vi« (vi,t) , (11) 



2V-1 N 



N! {Ni - 1) 



J2 E Wvfc-Vil, (8) 



fc=l j=k+l 



and Ao = v c tq, the characteristics particle mean- free- 
path. 



A. Thermodynamical limit 



being K (v vi, v 2 ) the probability that a given particle 
of the collision pair with initial velocities (vi, v 2 ) to reach 
v as a final velocity: 



Jf(v| v llV2 )= r ^[v-v(vi,v 3i fl)], (12) 
■/o o- 

being 



v(vi,v 2 ; 



vi + v 2 



+ 



vi - v 2 



n(# o + 0), (13) 



It could be spoken about a continuous distribution of 
the particles velocities when N tends to infinity because 
of the particles fill densely all the accessible velocities. 

Let n(v,r)<i 2 v be the relative density of those parti- and do (6) the differential cross section for the binary 



clcs remaining in the system whose velocities are located 
inside the differential region (v, v + dv) in a given com- 
putational step t. The function n (v, r) represents the 
relative density distribution function of particles in the 
neighborhood of the point v in the space of admissi- 
ble velocities E = {v £ R 2 : v 2 < i; 2 }. The variation 



encounters. The difference between Ii n and I ou t extracts 
out all those cases in which the interacting particles keep 
the same velocities after the interaction. 

The specific form of the irreversible dynamical equa- 
tion © leads to some straightforward consequences. Be- 
cause the probability p (vi , v 2 ) that a given particle of 
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the collision pair (vi,V2) to remain in the system after 
the collision: 



p(vi,v 2 ) = Jd 2 vK(v\ vi,v 2 ) : 



(14) 



may differ from the unity, < p(vi,v 2 ) < 1, the total 
relative density of the system p (r) = J E d 2 v n (v, r) will 
decrease or remain constant: 



— p(r) = -vjj d 2 v l d 2 v 2 q (vi,v 2 ) n (vi,r) n (v 2 ,r) < 0, 

SxE 

(15) 

being q (vi, v 2 ) = 1— p (vi, v 2 ), the escaping probability. 

A geometrical visualization of the inequality (|15|) is 
shown in FIG 121 Here X is the admissible subset of ve- 
locities in R 2 ; V represents the velocity of the center 
of mass, V = \ (vi + v 2 ), while u is the final velocity 
of the specific particle in the center of mass frame, that 
is, u = ~ |vi — v 2 | n (8q + 9). Only those values of the 
vector addition v = V + u belonging to £ correspond 
to the final state where the specific particle of the pair 
remains in the system after the collision. In the par- 
ticular case of the uniform differential cross section JSJ, 
p(vi,v 2 ) is numerically equal to fraction of the circum- 
ference perimeter described by the vector u belonging to 
the subset X. A simple geometrical analysis yields the 
following expression for this probability: 







when x < — 1 



p(vi,v 2 ) = < \ + i arcsin (x) when \x\ < 1 >, (16) 



when x > 1 



where x is given by: 



2|V||u| 



(17) 



It is evident that only those collision pairs satisfying 
the condition: 



1 

— Vl 

2 1 



i|vi-v 2 |<«c, (18) 



will remain in the system after the interaction with a 100 
% of probability, and this condition does not depend on 
the differential cross section used to describe the binary 
collisions, nor on the dimension of the space. 

As already shown, the relative density distribution 
function n (v; r) obeys a nonlinear irreversible dynamics, 
which is also dissipative. Thus, the evaporation is explic- 
itly taken into account in the dynamical equation of the 
system. The nonlinear nature of equation |J3J is a fea- 
ture making the difference from the standard treatment 
of evaporation by using a linear Fokker-Planck equation 
applied to the astrophysical systems 0, 0, |2(], |^| . This 
equation does not necessarily lead to exponential decays 
of the particles number with the time. 



B. Numerical experiments 

We set without lost of generality m = 1, Ao = 1 and 
u c = 2.6 in our numerical study of dynamics. The num- 
ber M of interacting pair chosen in a given computational 
step was settled to be M = [fiN], where fi = 0.05, in or- 
der to ensure a smooth temporal evolution of the phys- 
ical observables. We start from an initial configuration 
where all particles belong to the system, that is, bk = 1, 
Vfc £ [1,2, ...N]. Preliminary runs by using the uniform 
differential cross section (JSJ) allows us to set N — 1000, 
since size effects are not significant when we compare the 
results obtained by considering N — 500 and 1000. 

We are interested in describing some kind of dynamical 
evolution of the model system which could be considered 
as a quasi-stationary regime. To this aim we take into 
account the following initial distribution functions of par- 
ticle energies: 



Uniform Distribution (UNIC): 



m (e) = 



1 with < e < Ei < u c 
otherwise 



(19) 



• Truncating Boltzmann Distribution (TBIC 

/3exp( — (3e) 



mi (e) 



1 - exp (-/3u c ) 



1 with < e < u c 
otherwise 



(20) 



Double Size Uniform Distribution (DSIC): 



mil (e) 



Ei — Ei 



1 with Ei < e < Ei < u c 
otherwise 



(21) 

where E\, E% and j3 are positive real parameters. The 
particles velocities were uniformly distributed at random 
in all possible directions (0, 27r) for the three initial con- 
ditions. 

We start from an initial condition UNIC with Ei = 
2.0 by using the uniform differential cross section. FIG|3] 
shows the dynamical evolution of p and e in terms of 
the physical time t. It is easily seen that the system 
experiences a sudden relaxation process in which the 25 
% of the particles is evaporated in the first stages of the 
evolution. There, the system enters in a long transient 
regime where the evaporation events have been extremely 
reduced, but they still take place causing a slow variation 
of the kinetic energy per particle. 

The insert graph evidences that this relaxation process 
is progressive, that is, no real equilibrium configuration 
seems to be reached by extending the computational sim- 
ulation of the dynamics, but the rate of the evolution 
seems to be continuously reduced. 

Let us analyze now the effect of the particle evapora- 
tion in the energy distribution function of the system. 
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FIG. 3: Dynamical evolution of the relative density p and the 
kinetic energy per particle e by using the uniform differential 
cross section with UNIC initial condition with E\ = 2.0. The 
same dependences are shown in the insert graph by using a 
logarithmic scale for the time variable. 
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FIG. 4: Energy distribution function along the dynamical 
evolution shown in FIG[3] The system seems to stabilize in a 
Truncating Boltzmann distribution. 



We will apply the following procedure in order to ob- 
tain an appreciable statistics of the microscopic configu- 
rations around a given computational step r: We carry 
out certain number v of virtual trajectories, with 100 
computational steps each, by starting from the initial 
configuration at the computational step r. The data is 
collected by considering 20 microscopic configurations in 
each of such virtual trajectories, and therefore, a total 
of 20t> microscopic configurations are taken into account 
in order to build the energy distribution function. This 
procedure is applied in the successive analysis. 

FIG 0] shows the energy distribution function by con- 
sidering 100 virtual trajectories: at the first stages of the 
evolution with r = 0, 500, 1000, and a final configuration 
when t — 13000. It is remarkable that after a very short 
relaxation, the energy distribution function is stabilized 
in a truncating Boltzmann distribution of the form given 
in the equation (|2(}|> . which is easily evidenced in the 
linear dependence of the distribution functions shown in 
this figure by using a logarithmic scale. The existence of 
this profile allows us to claim that the system exhibits 
a quasi- stationary evolution where the relative density 
p and the energy per particle e seem to be the macro- 
scopic control variables along the successive dynamics of 
the system. 

Let us now perform a comparative study in order to 
clarify the possible influence of the initial conditions in 
the nature of the quasi-stationary evolution of the sys- 
tem. We consider three initial conditions for this analy- 
sis: a) UNIC: E x = 2.0; b) TBIC: = 0.25; c) DSIC: 
Ei = 1.0 and E 2 = 2.0. The results are shown in FIGE3 



and FIG HO As expected, the particle evaporation de- 
pends strongly on the initial microscopic configurations. 
However, the energy per particle e quickly tends to a 
regime where the successive evolution is so slow that the 
dynamical evolutions of this observable starting from dif- 
ferent initial microscopic configurations converge asymp- 
totically all them. 

FIG [7] and FIG|H1 show the energy distribution func- 
tions for these initial conditions at two different compu- 
tational times by considering also 100 virtual trajectories: 
FIGdat r = 1000, and FIG|H|at r = 13000. These re- 
sults strongly suggest us that the Truncating Boltzmann 
distribution H20|) is the quasi-stationary profile appearing 
in the dynamical evolution of a gas of particles driven by 
binary encounters which undergoes particle evaporation 
by using the uniform differential cross section. The great 
similarity among these profiles is a consequence of the 
sudden evolution of the system towards configurations 
whose energy per particle is around e ~ 0.27. For this 
energy the evaporation rate is so low that the subsequent 
dynamics spends too much time for exhibiting a signifi- 
cant variation. 

All those interacting particles whose kinetic energies 
are less than u c /2 do not provoke evaporation events, 
and therefore evaporation events could be only produced 
when a particle, whose energy e belongs to the energetic 
range ^u c < e < u c , is involved in a binary encounter. 
This is precisely the origin of the slow dynamical evolu- 
tion observed in the above numerical experiment when 
e ~ 0.27: Most of the particles remaining in the system 
belong to the stable region < e < \u c - 

As already noted, the uniform differential cross section 
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FIG. 5: Dynamical evolution of the relative density by using 
three different initial conditions and the uniform differential 
cross section. 



FIG. 7: Energy distribution function at r = 1000 for the 
dynamical evolution shown in FIGQ>]and FIG|S| 
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FIG. 6: Dynamical evolution of the energy per particle by 
using three different initial conditions and the uniform differ- 
ential cross section. 
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FIG. 8: The same observable shown in FIG0 but this time 
at t = 13000. 



(J5J allows a given interacting pair to explore; in an ef- 
fective fashion, all those possible results of their collision 
event, and therefore, large variations of the kinetic ener- 
gies of the interacting particles are very likely. We are 
interested in analyzing if the quasi-stationary profile ob- 
tained above remains unchanged even when the collision 



events are not so effective, that is, when large variations 
of the kinetic energy of the collision particles are not so 
favored. This analysis can be performed by using the 
differential cross section Q with a = 0.2ir. As already 
shown in FIG^ this particular value produces a moder- 
ate asymmetry among large and short deflections angles 
during the binary encounters, and therefore, a weaker ef- 
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FIG. 9: Dynamical evolution of the relative density p when 
a moderate asymmetry among large and short deflections is 
considered. 
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FIG. 11: Energy distribution function at r = 1000 for the 
dynamical evolutions shown in FIG [5] and FIG I10I 
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FIG. 10: Dynamical evolution of the energy per particle e by 
using a moderate asymmetry among large and short deflec- 
tions. 
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FIG. 12: The same observable shown in FIG I11I but this time 
at t = 13000. 



fectiveness of the equilibration mechanism with regard to 
the uniform case. 

We implement the same initial conditions used in the 
previous study. The dynamical evolution of the relative 
density p and the energy per particle e are shown in FIG|5] 
and FIG I1UI respectively. Despite we have considered a 



new microscopic interacting picture for the binary en- 
counters, it is evident here that no significant qualitative 
differences are exhibited in the dynamical behavior of 
these observables. This conclusion is reinforced by look- 
ing at the results showed in FIGI11I and FIG^l where 
it can be appreciated the persistence of the truncating 
Boltzmann profile l|2l)|l even after using a moderate asym- 
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FIG. 13: Dynamical evolution of the relative density p by 
using a differential cross section with a great asymmetry be- 
tween short and large deflections. 
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FIG. 15: Energy distribution function at r = 1000 for the 
dynamical evolutions shown in FIG I13l and FIG I14I 
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FIG. 14: Dynamical evolution of the energy per particle e FIG. 16: The same observable shown in FIGlHl but this time 
by using a differential cross section with a great asymmetry at t = 13000. 
between short and large deflections. 



metry among large and short deflections, although there 
is a slight tendency to a reduction of the particles pop- 
ulation for the large energies with a significant incidence 
of fluctuations. 

Let us now consider a very asymmetric differential 
cross section. We repeat the same numerical experiment 



described above, but this time by considering a = 0.027T. 
This value of the deformation parameter a leads to a pre- 
dominance of the small variation of the kinetic energies 
of the collision particles, and therefore, a significant re- 
duction of the effectiveness of the microscopic dynamics 
in reaching their quasi-equilibrium conditions. 

The evolution of p and e are shown in FIG I13I and 
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FIG. 17: Energy Distribution function with a = 0.02-7T at 
different energy per particles with their fitting by using the 
Michie-King like profile. The insert plot shows the energy 
dependence of the parameter ec- 



FIG. 18: Energy Distribution Function when e = 1.0 for dif- 
ferent values of deformation parameter a. The insert plot 
shows the dependence ec versus a. 



FIG I14I In a similar way, the energy distribution func- 
tions when t = 1000 and 13000, in FIG[]31and FIGlTol 

is also shown after running 200 virtual trajectories. As 
expected, the rate of evaporation has been reduced. The 
energy distribution function clearly evidences a new qual- 
itative effect of the evaporation, which is more accentu- 
ated in FIG El The results shown in these figures evi- 
dence a significant deviation of the distribution function 
from the linear dependence close of the energy cutoff u c , 
without vanishing at this energy. 

Let us now investigate, qualitatively, these quasi- 
stationary profiles by using a Michie-King-like distribu- 
tion : 

UJ M -K (e) = C (exp [-/3e] - exp {-pec)), (22) 

with < e < u c , being u c < ec < +oo. It is easy to 
note that when ec decreases from the infinity towards 
the cutoff energy u c , the energy parameter ec causes 
a continuous deformation of the Truncating Boltzmann 
profile l|20[) towards the well-known Michie-King profile 
of the globular clusters [T^ . l2lj . 

We perform several runs and obtain the energy dis- 
tribution functions at different energies, starting from a 
TBIC initial condition with (3 = 0.01. These results are 
shown in FIG 1 171 The reader can notice that the profile 
(|22[) offers a fairly good fitting for all these dependences, 
although for the case where e = 0.4 some discrepancies 
are evident. Such differences probably originate from the 
method used to obtain the distribution functions. The 
insert plot suggests a weak dependence of the parameter 
ec on the energy per particle e. 



We repeat the calculations by considering different val- 
ues of the deformation parameter a and keeping the en- 
ergy per particle fixed at e = 1.0. These results are shown 
in FIG^O It is clearly evidenced how the energy pop- 
ulation at large energies is increased when a decreases, 
which could be quantitatively described by the decreas- 
ing of the parameter ec- The origin of this behavior 
follows from the fact that the evaporation events take 
place practically when the involved particles are those of 
kinetic energy close to the scape energy u c . As already 
noted, the effect of the correction to the exponential pro- 
file is more noticeable with the increasing of the energy 
per particle e. Apparently, the specific values of e for the 
profiles shown in FIGlllland FIG ^1 do not allow us to 
observe this effect, which is more evident for the results 
shown in FIGlTgl 

Summarizing: The above numerical experiments allow 
us to conclude that the gas driven by conservative binary 
encounters and particle evaporation follows certain quasi- 
stationary evolution which seems to be macroscopically 
controlled by the energy per particle and the particles 
density. The generic form of the quasi-stationary profiles 
is practically independent of the initial conditions, be- 
ing appropriately describe by a Boltzmann profile trun- 
cated at the escape energy when the binary encounters 
favor an effective exploration of all velocities stable un- 
der the evaporation. The reduction of the effectiveness 
of this equilibration mechanism produces a continuous 
deformation of the truncating Boltzmann profile towards 
the well-known Michie-King distribution |19l |21| , which 
could be qualitatively described by the profile shown by 
equation 1221) . 
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III. GENERIC TREATMENT 

Let us consider the evaporation of the gas of binary 
encounters in a more general framework. Let P m be 
the probability that a given particle be in the state m, 
where m represents certain admissible velocity. The rate 
of change of this probability, JC m (P): 

^ = JC m (P) , (23) 
ar 

can be divided in two terms as follows: 

K m (P) - (Eq.Trans) m + (NEq.Trans) m . (24) 

The term (Eq.Trans) m considers all those binary encoun- 
ters affecting the state m without provoking evaporation 
events: 




0.5 1 .O 1.5 2.0 



energy 



(Eq.Trans) m = ^ w (AB)( m fa) (PaPb - P m Pfh) , 



ABr. 



(25) 

while (Evap. Trans) m involves all those encounters where 
a collision particle escapes from the system: 



(Evap.Trans) m = ^ D( 

(sm)t-Ps-Prn • 

ij st 

(26) 

It is very easy to note that the transition probabil- 
ity W / ( J 4s)(mm) °f t ne first term involves four admissible 
states, while the transition probability D^ m of second 
one only involves three states. This is the origin of the 
particles and energy dissipation by evaporation events, 
which are described by the equations: 



ar ar 



(27) 



being P = J2 m p ™ and E = J2 m £ mP m , where e m repre- 
sents the energy of the particle at the m-th state, as well 
as: 



£(P) 



^ t D(ij)mPiPj > 



and 



TZ (P) = - J2 DmmPiPj (e< + £j - Em) 



(28) 



(29) 



Let us now introduce the normalized distribution func- 
tion f m = P m /P ( E m /m = 1 ), which obeys the dy- 
namical equations: 



df m 
dr 



P [Km (/) " /m/C (/)] 



(30) 



FIG. 19: The plane s versus e reveals us that the model sys- 
tem follows certain sequence A-B along its dynamical evo- 
lution which does not depend on the initial conditions, be- 
ing this a clear manifestation of its quasi-stationary charac- 
ter. The curve A'-B' is the quasi-stationary sequence when 
a — 10.07T. The thick arrow represents the direction of the 
dynamical evolution along these sequences. 



where !C m (/) and /C (/) were defined at the equations 
(EU and 

Two important macroscopic quantities of this system 
are the energy per particle e [/] and the entropy func- 
tional s [/]: 



[/] = ^2 £ ™f™, S [/] = - ^ fm In fn 



(31) 



Since the distribution function f m must be microscop- 
ically controlled by the energy per particle e during the 
quasi-stationary regime, the entropy s [/] versus e plane 
must exhibits certain sequence revealing a functional de- 
pendence among these macroscopic quantities. FIGlT^l 
confirms this supposition by using three different initial 
conditions and considering a = 0.02-7T in the differential 
cross section (@J. 

The reader can notice that after a fast relaxation the 
entropy s evolves along certain curve A-B in the plane s 
versus e which does not depend on the initial conditions. 
This is a clear manifestation of the quasi-stationary char- 
acter of the dynamical evolution of the model system. 
Despite the reduction of the effectiveness of the equili- 
bration mechanisms for low values of the deformation 
parameter a, the quasi-stationary character of the dy- 
namical evolution of the model system seems to very ro- 
bust. For comparison, we also show in this figure the 
quasi-stationary sequence when a = 10.07T. It is easy to 
see that the decreasing of the deformation parameter a 
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provokes an increasing of disorder for a given energy per 
particle, which is directly related with the increasing of 
the particles population at large energies (see FIG lT%|l . 

The energy per particle e [/] evolves according to the 
equation: 

±e[f} = P{-K[f}e[f]+K[f}}, (32) 

being 1Z [/] the functional of the equation (|29[1 . while the 
dynamics of the entropy functional s [/] is described by 
the equation: 

^-s [/] = P {-JC [/] e [f] + W [/] +V[f}}, (33) 
where W [/] is given by: 

W Lfl = \ W (AB)(mfh) [In (/a/s) - In (fmfrh)} 

ABmfn 

X (f A f B - fmfrh) , (34) 

and V [/] as: 

V [/] = ^ D (ij) mfifj (In /i + In - In / m ) . (35) 

It is easy to see that the term W [/] > 0, where the 
identity only takes place with the establishment of the 
detailed balance: 

/(vi)/(v 2 ) = /(v 1 )/(v 2 ), (36) 

among all those pairs of velocities (vi,V2) and (vi,V2) 
belonging to the admissible space £ which are related by 
a collision event without evaporation. On the other hand, 
the functional V [f] does not exhibit, generally speaking, 
a definite signature. The energies e mi and Sj corre- 
sponding to the states related by an evaporation event, 
(ij — > m) are ordered as follows: e m < min(ei,£j). It is 
not difficult to show that T> [/] < (or V [/] > 0) when- 
ever / (e) be a decreasing (or increasing) monotonic func- 
tion on e, requirement satisfied by the Michie-King-like 
profile 

The quasi-stationary regime is established by the com- 
petition of two tendencies: First: the evolution towards 
the most likely macroscopic configuration where the de- 
tailed balance is imposed, which is represented by the 
presence of the conservative term W [/] in the equation 
(|33|l : Second one: the tendency of finding more stability 
under the evaporation, related with the presence of the 
dissipative terms — K. [/] s [f] + T>[f]. The functionals 
K> [/L ^ [/] an d Tl [/], despite being independent, vanish 
simultaneously when a monotonic distribution function 
on the energy / (e) , which drops to zero in the range 



u c /2 < e < u c , is considered. Since the stability of these 
tendencies cannot be simultaneously satisfied, the system 
evolves towards certain intermediate distribution which 
is never stationary. 

It is well-known the difficulties to find an appropriate 
criterion leading to the theoretical prediction of quasi- 
stationary profiles. We only found an appropriate crite- 
rion when the system equilibration mechanism is effective 
in exploring all accessible configurations which are stable 
under the evaporation: the system practically reach the 
most probable one, being this configuration the quasi- 
stationary profile. In those cases, the contribution of 
the dissipative terms in the entropy dynamical evolution 
seems to be only comparable to the term W [/] for con- 
figurations close to the most likely because of the pre- 
dominance of binary transitions which does not involve 
particle evaporation. However, the picture turns out to 
be considerably difficult in the general case with the re- 
duction of the effectiveness of the equilibration mecha- 
nisms. 

As already stressed, during the dynamical evolution 
the system distribution functions will converge towards 
certain sequence which is identified with the quasi- 
stationary one. Once settled at a given point of this 
sequence, the system will evolve along it. If we know 
some configuration belonging to this sequence, the pro- 
gressive evolution of the system will reproduce the rest 
of the quasi-stationary profiles. It is not difficult to un- 
derstand that this quasi-stationary sequence must con- 
tain the microscopic configuration with the maximum 
energy per particle e, f m (e). Since < e < u c , the 
distribution function with maxima energy per particle is 
fm (e) = 8 (e — u c ). Thus, in our opinion, the problem of 
finding the quasi-stationary profiles can be reduced to the 
determination of the temporal dependence of the distri- 
bution function / (e; r) by using the dynamical equation 
(|3l]|l with / (e; 0) = 5 (e — u c ), but this problem does not 
have a general analytical solution because of the nonlin- 
ear nature of the dynamics. 

We use the above dynamical criterion to build what 
we consider as the quasi-stationary sequence s qs versus 
e, which is shown in FIGlSOlby considering a = 0.027T. We 
obtain the data from two dynamical experiments, as well 
as by the extrapolation to the low energies region. The 
first experiment was performed by considering fi — 0.005 
in order to reduce the rhythm of the dynamical evolution 
when e is close to u c , and the second one, by considering 
fi = 0.1. The reader can recall that this parameter con- 
trols the temporal increase for each computational step, 
but it does not affect the dependence of distribution func- 
tion on the energy per particle. 

Both experiments were performed by considering N = 
20000 and by dividing the range < e < u c in 400 in- 
tervals in order to build the distribution functions. The 
function f m (e) was approximated by using a DSIC with 
E2 = uc and E± — u c — u c /400. The second experiment 
starts from a TBIC with f3 = —10.0, which exhibits a 
fast convergence towards the quasi-stationary sequence 
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FIG. 20: 



energy 



versus e dependence for the quasi-stationary 



evolution obtained by using the dynamical criteria when a = 
0.027T. The point A corresponds to the starting configuration 
for the second experiment. The insert plot shows the linear 
dependence at low energies for the s versus In e. 



of the first experiment. Notice that in this figure there 
is a superposition of both experiments in the interme- 
diate energetic range and no discrepancy is evidenced. 
The range < e < 0.22 was obtained from the extrapo- 
lation of the experimental data because of the evolution 
is extremely slow in this interval, and basically, the dis- 
tribution functions do not differ significantly from the 
truncating Boltzmann profile, yielding s ~ 6.035 + In e in 
this energetic range. 



IV. NONHOMOGENEOUS CHARACTER 
THE CONFINING POTENTIAL 



OF 



The model system considered in the above sections dis- 
regards the effect of the nonhomogeneous character of the 
confining potential energy. We are interested in describ- 
ing this effect, at least, in the extreme case where the 
system equilibration mechanism is very effective. As al- 
ready observed in the numerical experiments, when the 
particles are able to explore in an effective fashion all 
those microscopic configurations which are stable under 
the evaporation, the progressive evolution of the system 
leads to the setting of the ordinary equilibrium condi- 
tions, as an example, the imposition of the detailed bal- 
ance in the framework of the binary encounters. 

In the context of Hamiltonian systems, the equilibra- 
tion mechanisms relay on the strong chaoticity of the mi- 
croscopic dynamics. This property usually leads to the 
ergodicity and the equilibrium conditions are ordinarily 
described by considering a microcanonical distribution 



[rH| . If evaporation is present, the above picture can 
be naturally generalized by considering a microcanonical 
distribution where it has been disregarded all those mi- 
croscopic configurations where the particles are able to 
escape. The above argument was used in a very recent 
paper |22| in the context of astrophysical systems in or- 
der to develop an alternative version of the isothermal 
model of Antonov pj . 

The stable configurations for a Hamiltonian model of 
the form: 



H N =K N + V N 



N 



k=l 



Pi 



N-l N 



i—1 j>i 



n), (37) 



is the subset of the N-body phase space where the parti- 
cles satisfy the inequalities: pj./2 + u (r^) < £5, being es 
the energy cutoff, p/c and u(rk), the linear momentum 
and potential energy for the k-th particle: 



u(r k ) 



(38) 



where k = 1, ...N. We are interested in the specific form 
of the quasi-stationary distribution function / (r, p) in 
such conditions. 

This question was solved in ref.[23| in the context of 
astrophysical systems by considering the mean field ap- 
proximation appearing when N tends to infinity. The 
generalization of these results is straightforwardly fol- 
lowed by taking into account the generalized form of the 
mean field potential energy u (r): 



"M=/^(r-*)p(x). 



(39) 



where p (r) is the particles density at the point r of the 
physical space. Therefore, we will only expose here the 
final results. The interested reader can see the ref.[22] for 
details of this analysis in the framework of the Newtonian 
interaction. 

It can be proved that the particles density p (r) is given 
by the expression: 



/ m \ 2 r 

P (r) = (2^) ex P(-M + C (r) - (3u (r))F [y/¥(rj 

(40) 

where (3 and p are the canonical parameters which are 
specified by the energy and particles number constraints, 
F [z] is defined by the integral: 



F [z] = ir a / 4irx 2 exp (~x 2 ) dx, (41) 
Jo 

and $ (r) = (3 [es — u (r)]. The function C (r) is obtained 
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in a self- consistent way as follows: 



C (r) = / d d x 4> (r - x) p (x) 



<9u (x) 



•InF 



(42) 

Since p (r) = J <i 3 p / (r, p), it is very easy to note that 
equation (|40|) is derived from the distribution function 
/(r,p): 



/(r,p)=C exp(C [r)-/3e(r,p)] 



(43) 



which vanishes when e(r,p) > eg, being e (r, p) = 
^■p 2 + u (r). This expression differs from the truncating 
isothermal profile only because of the presence of the self- 
consistent function C (r) , which appears as an additional 
consequence of the particle evaporation due to the nonho- 
mogeneous character of the global interaction. Although 
the p versus u (r) dependence 140JI unifies the isothermal 
with the polytropic dependencies: p ~ exp (C + <&) when 
$ is large enough, while p ~ exp (C) <I>2 when $ << 1, 
the presence of the function C introduces a sensitive mod- 
ification to the features of the solutions. 

FIGl21lshows the effect of the function C (r) in the par- 
ticles density for the Newtonian potential 4>{r) = —n/r 
for configurations with spherical symmetry (details of 
this study are also found in ref.j^). The profiles A and 
C were obtained by considering the distribution function 
(|43|1 : profile A corresponds to a quasi-stationary configu- 
ration with low energy, which is characterized by the ex- 
istence of an isothermal core and a polytropic halo, while 
C corresponds to a high energy configuration where the 
isothermal core is not present. Profiles B and D corre- 
spond respectively to an isothermal and polytropic con- 



figurations with polytropic index 7 = 



polj 

E3 



(character- 



ized by the dependence p ($) oc $2). The presence of 
the function C (r) leads to a significant concentration of 
the particles toward the inner regions which goes beyond 
what the isothermal or the polytropic models predict. 

It is easy to understand that the function C depends 
in a self-consistent way on the energy function e(r, p). 
Therefore, the distribution function / (r, p), as a whole, 
depends only on the function e(r,p), that is, /(r,p) = 
J r [e(r, p)]. Needless to say that this is a requirement 
of the well-known Jeans theorem [17|. However, the 
consideration of such effect only starting from this last 
theorem is very unlikely. The present analysis suggests 
that the simple generalization of the homogeneous dis- 
tribution functions / (r, p) = ip (^p 2 ), by substituting 
2mP 2 ~~ * 2^P 2 +m (r) may be extremely naive, and prob- 
ably, most of these simple extensions are disregarding im- 
portant effects of the nonhomogeneous character of the 
global interaction in presence of particle evaporation. 

Thus, the consideration of a weaker chaoticity of the 
microscopic dynamics and the nonhomogeneous charac- 
ter of the interactions involve serious difficulties in the 
theoretical prediction of the quasi-stationary profiles. 
The reader may notice that the approach developed in 




1 E-3 



FIG. 21: The effect of the function C (r) in the particles den- 
sity for the Newtonian potential. 



ref.[23| and the standard treatment of some astrophys- 
ical structures, such as globular clusters and elliptical 
galaxies, by using the Michie-King profile [T^ . l2l| . are 
extreme alternatives to consider the effect of the evap- 
oration in such a framework. The applicability of such 
models depends crucially on how effective is the equili- 
bration mechanism in an actual astrophysical systems. 

It is well-established that the picture of the binary en- 
counters provide a nice justification for the use of the 
King's models. The basic explanation is found in the 
fact that the long-range character of the gravitational in- 
teraction favors much more the small energy interchange 
among those distant particles in comparison to the large 
ener gy interchanges occurring during the close encoun- 
ters iii m m . 

Nevertheless, the binary encounter is not the only equi- 
libration mechanism appearing in this context. In fact, 
it has been also shown the existence of a very strong 
chaoticity in gravitational N-body systems [Tj, which 
acts in time scales much smaller than the characteristic 
time for the close binary collisions, as well as it oper- 
ates for every bound energy. In the viewpoint of the au- 
thor of such work, this result not only confirms that "... 
bound collisionless self gravitating systems are character- 
ized by a dynamical instability proceeding at a very fast 
rate, but also suggests that they possess the strongest sta- 
tistical properties, in analogy with those of standard dy- 
namical systems in the regime of fully developed stochas- 
ticity ..." As expected, this argument supports the 
consideration of the ergodic hypothesis in this context, 
and therefore, the consideration of a suitably regularized 
microcanonical description for actual astrophysical sys- 
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terns USE 03. 

Despite the differences concerning to the specific form 
of the above energy profiles, the question is that the fea- 
tures of their corresponding spatial particle distributions 
possess more similarities than differences: both distribu- 
tions lead to finite configurations with isothermal cores 
and polytropic haloes, even purely polytropic configura- 
tions are also possible, where the system size is deter- 
mined from the tidals interactions. 

In the remarkable paper |l9j |. King performed a nu- 
merical comparison of the Michic-King model with the 
isothermal distribution truncated at the escape energy, 
without the consideration of self-consistent function C, 
which is shown in FIG. 2 and FIG. 7 of this reference. 
While the King model provides an excellent fit for the 
data corresponding to a high concentrate globular clus- 
ter, the truncated distribution is unable to fit an obser- 
vational data where a larger concentration of particles 
appears towards the inner region. In our view, this could 
be the very effect introduced by taking into consideration 
the self-consistent function C, but this last observation 
deserves a further analysis. Moreover, it is difficult to dis- 
tinguish the origin of such effect since there exist other 
physical backgrounds modifying the form of the particles 
distribution in the space, such as the mass segregation 
[2^ | (see also ref.[24| for a recent development). 

V. CONCLUSIONS 

The numerical experiments carried out for the gas of 
binary encounters clearly show the robustness of the 
quasi-stationary character of the dynamical evolution 
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